# This repository provides the implementation of our recent work:
# A Cost-Effective and Scalable Machine Learning Approach for Quality Assessment of Fresh Maize Kernel Using NIR Spectroscopy
# Authors: Shi et al.
# This code is released for academic and non-commercial use only.

#这是Select算法的迭代大循环
from __future__ import absolute_import,division,print_function
import numpy as np
from sklearn import metrics
import tensorflow as tf
import pandas as pd
from tensorflow import keras
import os
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
from tensorflow.keras.models import load_model

# 设置字体
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# Rice
# a=[13]
# b=['relu']
# c=[2100]
# d=[49]
# e=['relu']
# f=[2800]


# Wheat
# a=[6]
# b=['swish']
# c=[3000]
# d=[45]
# e=['swish']
# f=[3000]

# Barley
a=[42]
b=['swish']
c=[3200]
d=[40]
e=['swish']
f=[2600]

# 稻米数据集
# filename = 'A_Train_test/Rice_Train_test.xlsx'
# Index = pd.read_excel(filename, sheet_name="30_182")
# 小麦数据集
# filename = 'A_Train_test/Wheat_Train_test.xlsx'
# Index = pd.read_excel(filename, sheet_name="30_83")
# 大麦蛋白质数据集
filename = 'A_Train_test/Barley_Train_test.xlsx'
Index = pd.read_excel(filename, sheet_name="30_128")

train_non_empty_count = Index.iloc[:, 0].dropna().shape[0]
train_index = Index.iloc[0:train_non_empty_count, 0].astype(int) -1
test_non_empty_count = Index.iloc[:, 1].dropna().shape[0]
test_index = Index.iloc[0:test_non_empty_count, 1].astype(int) -1

# 稻米数据集
# df1 = pd.read_excel('Rice_cm/Interpolated_Rice_1557.xlsx')   #1557波长点
# 小麦数据集
# df1 = pd.read_excel('Wheat_cm/Interpolated_Wheat_1557.xlsx')   #1557波长点
# 大麦蛋白质数据集
df1 = pd.read_excel('Barley_cm/Interpolated_Barley_1557.xlsx')   #1557波长点

df2 = pd.read_excel('Sythentic/Synthetic.xlsx')
X_R = df1.iloc[:, 1:].reset_index(drop=True)
Y_R = df1.iloc[:, 0].reset_index(drop=True)
X_S = df2.iloc[:, 1:].reset_index(drop=True)
Y_S = df2.iloc[:, 0].reset_index(drop=True)
Y_S = Y_S * 0.92    # 蛋白质纯度为0.92

X_S.columns = pd.to_numeric(X_S.columns, errors='coerce')
X_R.columns = pd.to_numeric(X_R.columns, errors='coerce')

X_R_train = X_R.iloc[train_index, :]
X_R_test = X_R.iloc[test_index, :]
Y_R_train = Y_R.iloc[train_index]
Y_R_test = Y_R.iloc[test_index]

# 处理方式1
X_all = pd.concat([X_S, X_R_train, X_R_test], axis=0).reset_index(drop=True)
Y_all = pd.concat([Y_S, Y_R_train, Y_R_test], axis=0).reset_index(drop=True)

num_components = 26
pca = PCA(n_components=num_components, svd_solver='full')
pca.fit(X_all)
mm = MinMaxScaler()
mm.fit(Y_all.values.reshape(-1, 1))

for k in range(len(a)):
    tf.random.set_seed(47)

    X_train1 = pca.transform(X_S)  # 第一个网络的训练集，人工数据
    Y_train1 = mm.transform(Y_S.values.reshape(-1, 1))
    X_test11 = pca.transform(X_R_train)
    Y_test11 = mm.transform(Y_R_train.values.reshape(-1, 1))
    X_test12 = pca.transform(X_R_test)
    Y_test12 = mm.transform(Y_R_test.values.reshape(-1, 1))

    model1 = keras.models.Sequential([
        keras.layers.Dense(a[k], activation=b[k],
                           input_shape=X_train1.shape[1:],
                           kernel_initializer=tf.random_normal_initializer(),
                           bias_initializer=tf.constant_initializer(0.0)),
        keras.layers.Dense(1,
                           kernel_initializer=tf.random_normal_initializer(),
                           bias_initializer=tf.constant_initializer(0.0)),
        ])
    model1.summary()
    model1.compile(loss='mse', optimizer = keras.optimizers.Adam(learning_rate=0.001))
    history1 = model1.fit(X_train1, Y_train1,
                       epochs =c[k])

    model1.save("model1.h5")
    my_model1 = load_model("model1.h5")

    Y_pre1 = my_model1.predict(X_train1)  #训练集预测结果
    Y_pre11 = my_model1.predict(X_test11)   #测试集预测结果
    Y_pre12 = my_model1.predict(X_test12)   #测试集预测结果

    #上面得到的是由第一个神经网络预测真实玉米数据出来的结果
    Y_pre1 = Y_pre1.reshape(-1)
    Y_pre11 = Y_pre11.reshape(-1)
    Y_pre12 = Y_pre12.reshape(-1)
    #****************************下面准备搭建第2个网络*****************
    result1 = []
    for i in range(len(Y_pre11)):
        result1.append(Y_pre11[i] -Y_test11[i])
    print("result1:",result1)
    result1 = [-x for x in result1]
    result2 = []
    for i in range(len(Y_pre12)):
        result2.append(Y_pre12[i] -Y_test12[i])
    print("result2:",result2)
    result2 = [-x for x in result2]

    #训练集
    X_train21 = X_test11
    Y_train21 = result1
    Y_train21 = np.array(Y_train21)
    Y_train21 = Y_train21.astype(np.float64)
    #测试集
    X_train22 = X_test12
    Y_train22 = result2
    Y_train22 = np.array(Y_train22)
    Y_train22 = Y_train22.astype(np.float64)

    model2 = keras.models.Sequential([
        keras.layers.Dense(d[k], activation=e[k],
                           input_shape=X_train21.shape[1:],
                           kernel_initializer=tf.random_normal_initializer(),
                           bias_initializer=tf.constant_initializer(0.0)),
        keras.layers.Dense(1,
                           kernel_initializer=tf.random_normal_initializer(),
                           bias_initializer=tf.constant_initializer(0.0)),
    ])
    model2.summary()
    model2.compile(loss='mse', optimizer = keras.optimizers.Adam(learning_rate=0.001))
    history2 = model2.fit(X_train21, Y_train21, epochs =f[k])
    model2.save("model2.h5")
    my_model2 = load_model("model2.h5")

    Y_pre21 = my_model2.predict(X_train21)  #训练集预测结果
    Y_pre22 = my_model2.predict(X_train22)  #训练集预测结果
    Y_pre21 = Y_pre21.reshape(-1)
    Y_pre22 = Y_pre22.reshape(-1)

    predict1 = []
    for i in range(len(Y_pre11)):
        predict1.append(Y_pre11[i] + Y_pre21[i])
    print("predict1:",predict1)

    predict2 = []
    for i in range(len(Y_pre12)):
        predict2.append(Y_pre12[i] + Y_pre22[i])
    print("predict2:",predict2)

    # 下面进行反归一化
    predict1 = np.array(predict1).reshape(-1, 1)
    predict2 = np.array(predict2).reshape(-1, 1)
    predict1 = mm.inverse_transform(predict1)
    predict2 = mm.inverse_transform(predict2)
    predict1 = np.array(predict1).reshape(-1)
    predict2 = np.array(predict2).reshape(-1)

    Y_train1 = np.array(Y_S).reshape(-1)
    Y_test11 = np.array(Y_R_train).reshape(-1)
    Y_test12 = np.array(Y_R_test).reshape(-1)

    #评价指标
    msec = metrics.mean_squared_error(Y_test11,predict1)
    rmsec = np.sqrt(msec)
    msep = metrics.mean_squared_error(Y_test12,predict2)
    rmsep = np.sqrt(msep)

    R1 = metrics.r2_score(Y_test11,predict1)
    R2 = metrics.r2_score(Y_test12,predict2)
    print('R1：', metrics.r2_score(Y_test11,predict1))
    print('R2：', metrics.r2_score(Y_test12,predict2))
    # 使用numpy库计算相关系数
    correlation1 = np.corrcoef(Y_test11,predict1)[0, 1]
    print('rmsec：',rmsec)
    print("Corr. Coeff.:", correlation1)
    correlation2 = np.corrcoef(Y_test12,predict2)[0, 1]
    print('rmsep：',rmsep)
    print("Corr. Coeff.:", correlation2)

    #下面计算RPD=S.D/RMSEP
    #print("总体标准差:", np.std(f_pred2))
    SD=np.std(Y_test12, ddof=1)
    #SD=np.std(f2_data)
    RPD=SD/rmsep
    #print("SD",SD)
    print("RPD:", RPD)

# # *******************以下为绘图部分2(fig1和fig2也反归一化***************************
    Y_pre1 = np.array(Y_pre1).reshape(-1, 1)
    Y_pre1 = mm.inverse_transform(Y_pre1)
    Y_pre1 = np.array(Y_pre1).reshape(-1)
    Y_pre11 = np.array(Y_pre11).reshape(-1, 1)
    Y_pre11 = mm.inverse_transform(Y_pre11)
    Y_pre11 = np.array(Y_pre11).reshape(-1)
    Y_pre12 = np.array(Y_pre12).reshape(-1, 1)
    Y_pre12 = mm.inverse_transform(Y_pre12)
    Y_pre12 = np.array(Y_pre12).reshape(-1)

    Y_train21 = np.array(Y_train21).reshape(-1, 1)
    Y_train21 = mm.inverse_transform(Y_train21)
    Y_train21 = np.array(Y_train21).reshape(-1)
    Y_train22 = np.array(Y_train22).reshape(-1, 1)
    Y_train22 = mm.inverse_transform(Y_train22)
    Y_train22 = np.array(Y_train22).reshape(-1)
    Y_pre21 = np.array(Y_pre21).reshape(-1, 1)
    Y_pre21 = mm.inverse_transform(Y_pre21)
    Y_pre21 = np.array(Y_pre21).reshape(-1)
    Y_pre22 = np.array(Y_pre22).reshape(-1, 1)
    Y_pre22 = mm.inverse_transform(Y_pre22)
    Y_pre22 = np.array(Y_pre22).reshape(-1)

    # fig1
    plt.figure(figsize=(6, 5))
    plt.scatter(Y_train1, Y_pre1, color='blue', label="Synthetic", s=40)
    plt.scatter(Y_test11, Y_pre11,  color='red', label="Calibration", s=40)
    plt.scatter(Y_test12, Y_pre12,  color='green', label="Validation",s=40)
    # 添加45度对角线
    min_val = min(np.min(Y_train1), np.min(Y_test11), np.min(Y_test12),np.min(Y_pre1),np.min(Y_pre11),np.min(Y_pre12))
    max_val = max(np.max(Y_train1), np.max(Y_test11), np.max(Y_test12),np.max(Y_pre1),np.max(Y_pre11),np.max(Y_pre12))

    plt.plot([min_val, max_val], [min_val, max_val], 'k--')  # 黑色虚线
    plt.xticks(fontsize=15)  # 设置 x 轴刻度标签的字体大小
    plt.yticks(fontsize=15)  # 设置 y 轴刻度标签的字体大小
    # 添加图例和标签
    plt.legend(loc='lower right', fontsize=15, frameon=False)  # 右下角无边框图例
    plt.xlabel('$y$(Reference)', fontsize=15)  # 使用 LaTeX 语法设置下标
    plt.ylabel('$y^{0}$', fontsize=15)  # 使用 LaTeX 语法设置下标
    # 在图中标注指标
    text_position_x = min_val + 0.01 * (max_val - min_val)
    text_position_y = max_val - 0.01 * (max_val - min_val)
    # 添加文本信息
    plt.text(text_position_x, text_position_y,
             # f'Rice(PCNN)\n',
             # f'Wheat(PCNN)\n',
             f'Barley(PCNN)\n',
             fontsize=15, verticalalignment='top', color='black')
    plt.tight_layout()
    # plt.show()

    # fig2
    plt.figure(figsize=(6, 5))
    plt.scatter(Y_test11, Y_pre21, color='red', label='Calibration', s=40)
    plt.scatter(Y_test12, Y_pre22, color='green', label='Validation', s=40)
    # 添加45度对角线
    min_val = min(np.min(Y_test11), np.min(Y_test12), np.min(Y_pre21), np.min(Y_pre22))
    max_val = max(np.max(Y_test11), np.max(Y_test12),np.max(Y_pre21), np.max(Y_pre22))
    plt.plot([min_val, max_val], [min_val, max_val], 'k--')  # 黑色虚线
    plt.xticks(fontsize=15)  # 设置 x 轴刻度标签的字体大小
    plt.yticks(fontsize=15)  # 设置 y 轴刻度标签的字体大小
    # 添加图例和标签
    plt.legend(loc='lower right', fontsize=15, frameon=False)  # 右下角无边框图例
    plt.xlabel('$y$(Reference)', fontsize=15)  # 使用 LaTeX 语法设置下标
    plt.ylabel(r'$y^{\epsilon}$', fontsize=15)  # 使用 LaTeX 语法设置下标
    # 在图中标注指标
    text_position_x = min_val + 0.01 * (max_val - min_val)
    text_position_y = max_val - 0.01 * (max_val - min_val)
    # 添加文本信息
    plt.text(text_position_x, text_position_y,
             # f'Rice(PCNN)\n',
             # f'Wheat(PCNN)\n',
             f'Barley(PCNN)\n',
             fontsize=15, verticalalignment='top', color='black')
    # plt.show()

    # # fig3
    # # 绘图
    plt.figure(figsize=(6, 5))
    plt.scatter(Y_test11, predict1, color='red', label='Calibration', s=40)  # 红色表示训练样本
    plt.scatter(Y_test12, predict2, color='green', label='Validation', s=40)  # 绿色表示测试样本
    # 添加45度对角线
    min_val = min(np.min(Y_test11), np.min(Y_test12), np.min(predict1), np.min(predict2))
    max_val = max(np.max(Y_test11), np.max(Y_test12), np.max(predict1), np.max(predict2))
    plt.plot([min_val, max_val], [min_val, max_val], 'k--')  # 黑色虚线
    plt.xticks(fontsize=15)  # 设置 x 轴刻度标签的字体大小
    plt.yticks(fontsize=15)  # 设置 y 轴刻度标签的字体大小
    # 添加图例和标签
    plt.legend(loc='lower right', fontsize=15, frameon=False)  # 右下角无边框图例
    plt.xlabel('$y$(Reference)', fontsize=15)  # 使用 LaTeX 语法设置下标
    plt.ylabel(r'${\hat{y}}$(Prediction)', fontsize=15)  # 使用 LaTeX 语法设置下标
    # 在图中标注指标
    text_position_x = min_val + 0.01 * (max_val - min_val)
    text_position_y = max_val - 0.01 * (max_val - min_val)
    # 添加文本信息
    plt.text(text_position_x, text_position_y,
             # f'Rice(PCNN)\n'
             # f'Wheat(PCNN)\n'
             f'Barley(PCNN)\n'
            f'$R^2_C$: {R1:.3f}    RMSEC: {rmsec:.3f}\n'
            f'$R^2_V$: {R2:.3f}    RMSEV: {rmsep:.3f}\n'
            f'RPD: {RPD:.3f}',
             fontsize=15, verticalalignment='top', color='black')
    plt.show()
    